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Abstract 

The self-energy screening correction is evaluated in a model in which the 

effect of the screening electron is represented as a first-order perturbation of 

the self energy by an effective potential. The effective potential is the Coulomb 

potential of the spherically averaged charge density of the screening electron. 

We evaluate the energy shift due to a ls^/a, 2s 1 / 2 , 2p]/ 2 , or 2p 3/ / 2 electron 

screening a ls]/ 2 , 2s]/ 2 , 2p 1 / 2 , or 2p 3 / 2 electron, for nuclear charge Z in the 

range 5 < Z < 92. A detailed comparison with other calculations is made. 
31.30.Jv 
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I. INTRODUCTION 



The self-energy correction to the electron-electron interaction is one of the many con- 
tributions of order o? to an atomic binding energy. These corrections, shown as Feynman 
diagrams in Fig. [I], are often called self-energy screening corrections, and for inner shells they 
are the largest of all fourth-order radiative corrections. They give rise to three terms, which 
are represented in Fig. ^ when one distinguishes between the reducible and irreducible part 
of the diagram in Fig. [I] (A). 

A first attempt to evaluate the contribution of such diagrams from bound state quantum 
electrodynamics (BSQED) was made in 1991 Jy] in an approximation in which the electrons 
not associated with the self-energy loop were represented as a perturbing potential (Fig. [|). 
The potential was obtained by taking a spherical average of the electron wave function, 
and calculating the potential associated with the resulting charge density. More recently, 
direct evaluations of the diagrams of Fig. |] for the ground state of two-electron ions have 
been made [0-0], and the method of H was used to provide the self-energy correction to 
several interactions ||. In this paper we report on a complete calculation of the self-energy 
screening in the approximation of |J, for all combinations of pairs of the states lsi/2, 2sy2, 
2pi/2> or 2^3/2 in the range 5 < Z < 92. 

Terms corresponding to the diagram in Fig. |3](A) are obtained by varying the self-energy 
expression with respect to the external potential of the bound-state wave function, leading 
to the introduction of the first-order correction to the wave function in the potential 5V. 
Terms corresponding to the diagram in Fig. ||](B) are obtained by varying the expression 
for the self-energy with respect to the external potential of the Green's function and while 
those in Fig. 0(A') correspond to variation with respect to the energy of the bound state. 

The expression for the self energy in a large class of potentials V(x) can be written as 
the sum £se = £l + £h of a low-energy part £l and a high-energy part £h given (in units in 
which h = c = m e = 1) by [[| 
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(£ n - Z) 2 X 2 1 
(1) 



and 

£k = tt^t / ete / dx 2 [ dxiipl l (x 2 )a a g(x 2l x 1 , z)a^ 1 ip n (x 1 )- 5m f dx <p n (x)/3tp n (x) 

Z7T1 JC K J J X 2 \ J 

(2) 

1 /2 

where b = — i [(£„ — z) 2 + i8\ , Re (6) > 0, and x 2 \ — x 2 — X\. In these expressions, tp n and 
S n are the eigenfunction and eigenvalue of the Dirac equation for the bound state n, and Q 
is the Green's function for the Dirac equation corresponding to the operator Q = (7i — z)^ 1 , 
where Ti = ct-p+V + f3 is the Dirac Hamiltonian. The indices I and m are summed from 1 to 
3, and the index /i is summed from to 3. The contour Ch extends from —zoo to — ie and 
from + ie to +ioo, with the appropriate branch of b chosen in each case. For the present 
calculation, we assume that the potential V(x) is close to a pure Coulomb potential, except 
for a small correction SV(x), which is not necessarily spherically symmetric. Indeed some 
applications of this method have been made with non-spherically-symmetric perturbations 
||. We obtain the screening correction to the self-energy by making the replacements 

V(x) -> V(x) + 6V(x) , (3a) 

ip n (x) -> <j) n (x) + 5(f) n (x) , (3b) 

G(x 2 , ail, z) G(x 2 , x 1 , z) + 8G(x 2 , x 1 , z) (3c) 

E n -> E n + 5E n (3d) 

in Eqs. (0) and (0) and retaining only the first-order correction terms. In the Eqs. (|l[) and 
(0), and (paf ) to (|3cJ) above, we use the symbols Q, S n , ip n for the exact quantities in the 
potential V(x), while the symbols G, E n , <fi n represent the corresponding exact quantities in 
a pure Coulomb potential V(x). The same conventions are employed throughout the paper. 
We denote the operations on the unperturbed self energy that lead to these three corrections 
by 5c/,, 5q, and 5e, respectively. In particular, we employ the notation 



d 

So = «?| (5) 
& = ^„A, (6) 



and the total correction is the sum 



d d d 

+ 5G^- + 5E n ^- (7) 



in which the partial differentiation symbol denotes formal differentiation with respect to the 
indicated variable, with the result evaluated with the unperturbed functions. 

In Sec. H we write expressions for the first-order perturbation corrections to the energy, 
the wave function and the Green's function. In Sees. |TI|, |V|, and we derive the expressions 
for the various contributions to the screened self-energy corresponding to the three diagrams 



of Fig. y. In a series of three earlier papers []TQ| [T2|| , we have derived and tested a method 



of analytically isolating divergent contributions to the self-energy diagram in coordinate 
space. In Sees. |V| and [V| we derive from this earlier work the generalizations of the analytic 
subtraction terms which are necessary to make all contributions to the self-energy screening 
finite. The numerical results are presented in Sec. [VT], and Sec. |V11| is the conclusion. 



II. PERTURBATION EXPANSION 

Replacing the potential V(cc) by V(x) + SV(x), where SV(x) is spherically symmetric, 
changes the wave functions, the energy, and the Green's function, which appear in Eqs. ([!]) 
and @. From standard perturbation theory the first-order energy correction is 

5E n = J dx<fil(x) 8V{x)(j) n (x) 

poo 2 

= dxX 2 J2fn,i( X )d V ( X )fn,i( X )> ( 8 ) 

J ° i=l 

where the radial wave function f n j is defined by writing 
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(9) 



and where Xk( x ) * s the Dirac angular momentum eigenfunction. 

The first-order correction to the wave function is given with the aid of the reduced 
Green's function G R (x 2 ,X\, E n ), defined by (see for example Ul3| ) 

4>m{x2)4L{xi) 



G R {x 2 , Xl ,E n )= J2 

m 



lim 

z^E n 



G(x 2 ,x 1 ,z) - 



E — E 

- LJ m J -'n 

<l>nM<f>U X l) 
E n - Z 



as 



^>n{x 2 ) = - J dXiG K {X2,X 1 ,E n )5V{Xi)(j) n {Xi) 



For a spherically symmetric potential we have 



(10) 



(11) 



dx x x\ G R>ij (x 2 , xi, E n ) SVix^fnjixt) 



(12) 



In Eq. ([12]) , the components of the radial reduced Green's function G R ' l ^(x 2 , xi, E n ) are 
defined in analogy with the components of the full Green's function G % i{x 2 , x\, z) as given 
in Ref. §, Eq. (A. 14). 

To evaluate the first-order correction to the Green's function we use the well-known 



expansion 



G{z) 



l 



H + SV-z 

1 6V- 1 



H — z H-z H-z 

5V^J— 6VrJ— + ■ ■ 



H-z H-z H-z 
and the term of first order in SV is 



(13) 



6G(z) 



1 sv- 1 



H-z H-z 



(14) 



which has second-order poles at the eigenvalues of the Dirac equation. In coordinate space, 
the first-order correction is 



SG(x 2 ,x 1 ,z) = 

- J dx 3 G(x 2 , x 3 , z) SV(x 3 )G(x 3 , x u z) , (15) 

and for SV(x 3 ) spherically symmetric, we have 

5G l Hx 2 ,xx,z) = — dx 3 x\ 
Jo 

xJ2Gi k (x 2 ,x 3 ,z)5V(x 3 )G%(x 3 ,x 1 ,z) . (16) 



k=l 



III. LOW-ENERGY PART 

The low-energy part, for an arbitrary external spherically symmetric potential, when 
integrated over the spherical angles of the vectors x 2 and Xi, yields 

S L = -S n + -P f £n dzU(z) (17) 

7T IT JO 

with 

oo 
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/*oo roo 

U(z) — — / dx 2 x\ \ dx\x\ 
Jo Jo 



2 

x Yj E ^n,i(x 2 )G t K 3 (x 2 ,x 1 ,z)J r nJ (x 1 )A^( y x 2 ,xi,z), (18) 

K i,j=l 

where the summation over k runs over all nonzero integers, and where i = 3 — % and J = 3 —j. 

We are concerned with the first-order perturbation in this expression that arises from 
variation of the external potential. The self energy depends on the potential through three 
quantities that appear in Eq. flTT| ) and (0), the wave function J-" n i (x), the energy eigenvalue 
£ n , and the Green's function Q l J(x 2 , x±, z). The three corrections are denoted by 8<f, £l, <5e £l, 
and <5g£l, respectively, with the total 

5E L = 5^£ L + 5 E £ L + 5 G £ L . (19) 
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A. Lower-order terms 



The expression ([!]) contains spurious parts of lower order in Za than the complete result. 
The physically significant part is isolated in a function jF L (Za) defined by ]I0| 



£1 



7T 



5 2 7 

+ - (Vnl^l^n) + T (Vn|V|v?n) 
D O 



(Za)' 



(20) 



Here the linear perturbation of this function with respect to variation of the external poten- 
tial is of interest. A corresponding function SF L (Za) is defined by 



6Et 



ft 



7T 



2 5E n + ^(<p n \P\S<Pn) + l(<Pn\V\5<j )n ) 



a(Za)'- 



■6F L (Za) 



where the fact that 



(21) 



5 (<p n \V\0 n ) = (<j> n \6V\<l> n ) + 2 (0„|V|50 n ) 
= 5E n + 2(<f> n \V\5<f> n ) 



(22) 



has been taken into account. 



B. Low-order matrix elements 

The lower-order expectation values involving the first-order correction to the wave func- 
tion S(p n can be evaluated by direct numerical integration. However, to get an independent 
check of the precision of the calculation, particularly when strong cancellation occurs, we 
derive a number of useful expressions in which 5<p n does not occur. The energy perturbation 
is given by the conventional expression 



5E n = (<f> n \6V\<j> n ) 







fiW+ffa) SV(x) . (23) 
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The wave-function perturbation terms are 



{(j) n \(3\5(j) n ) = dx 2 dxx4>\ l {x 2 )(3 



x 2^ — p _ p — (24) 



and 



(0 n |V|50 n ) = J dx 2 J dx 1 (/)l(x 2 )V(x 2 ) 



x 2^ F _ F — <*V(ari)0 n (xi) . (25) 

Since we are considering here the case where the unperturbed potential is the Coulomb po- 
tential with known wave functions, we can simplify the calculation of these matrix elements. 
In particular, we interpret the expressions in Eqs. and ( p5|) as perturbations of the wave 
function on the left-hand side to give 



») = (8 p <l> n \5V\<t> n ) (26) 
and 

where 5p<p n and 5y0 n are first-order corrections to the wave function due to perturbations (3 
and V, respectively. These are simply calculated as the coefficients of 5 in the power series 
expansions in 5 of the wave functions for the appropriately modified Hamiltonians 

Hf, = a-p+(l + 6)l3 + V 

= ct-p+{l + 5)/3m + V (28) 

and 

H v = a-p + (3+(l + 6)V . (29) 

In ( |2"8"D the second line restores the mass dependence of the Hamiltonian in order to exhibit 
the dependence of the modified Coulomb wave functions on 5. 
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The wave function correction 8p<f> n is obtained by replacing m by (1 + 8)m in n and 
calculating the coefficient of the term linear in 5, which is equivalent to writing 



9 A. 

m——<p n 
dm 



(30) 



and leads to 



(K\f3\S<p n ) 



m d 
2 dm 



(<Pn\5V\<l>n) 



(31) 



where the derivative acts only on the wave function. A further simplification is possible 



based on the mass dependence of the wave function in ([H]). If x is replaced by x/m in the 
second line of fl28|), then the mass factors out of the Hamiltonian, and the wave function is 
independent of m. As a result, we have 

m d r°° 



{4>nW<Pn) 



dx x 



2 dm Jo 



ft(x)+fi(x) SV(x) 



m 3 Jo 



x 



fi(x/m) + /| (x/m) 



m d 
2 dm 



— / dx x 3 
2 Jo 

-i(0 n |x(^)'|0 n } , 



SV(x/m) 
8 



(32) 



where the convention m = 1 has been restored in the last line, and (5V)' represents the 
operator corresponding to ^5V(x). 

The correction 5v<f> n is obtained by replacing Z by (1 + S)Z in (p n and calculating the 
coefficient of the term linear in 5 which is equivalent to writing 

„ d , 



>VVn 



dZ 



(33) 



which leads to 



(34) 



where it is understood that the derivative acts only on the wave function, or 



(<t>n\V\5<P n ) 



Z dx x 1 
Jo 



h{x) dJ ^ + H*) dh{x) 



dZ 



dZ 



SV{x) . 



(35) 



The derivative of the potential in (|3"2] ) is calculated analytically, and the derivatives in 
are calculated numerically by evaluating the wave function with 32 figure precision and using 
a symmetric derivative formula with 5(Za) = 10~ 12 . 

The matrix elements listed above are evaluated by Gaussian quadrature, and the code was 
tested in the Coulomb case and compared to the analytic results, as described in Appendix 

E 

C. Low-Energy energy-level and wave function correction 



The correction from the energy-level perturbation of the low-energy part of the self energy 
(0) ^ 



5 e £l = SE r , 



d£ r , 



(36) 



where 



a 


a 




TT 








-p/ 




TT JO 



dx 2 / dxi <f)l(x 2 )a l G(x 2 ,x 1 ,z)a m <p n (xi)(5i rn V 2 ■ Vi - V' 2 V" 



l^,^ S H( E n - Z)x 2 l] 



(E n - z) 2 X 21 



z=E„ 



E n 



The second term on the right-hand side of fl3"T|) makes no contribution because 



(E n - z) 2 x 2 i 3 



-5 lm (E n -z) + 0(E n - z) 



d sin[(E n - z)x 2 i 
dE n (E n - z) 2 x 2 i 



(38) 



and 



dx 2 [ dx x <Pi(x 2 )a l G(x 2 ,x 1 ,z)a m Mxi) = f dx 2 [ dx, <P n (x 2 )a l £ ^^^VUs.ji Oil) 



E n - z 



0(1). 



(39) 



10 



The first term on the right-hand side of Eq. ([39]) vanishes by virtue of the identity 

J dx (j>l(x)a l (j)j(x) = i J dx (j)l l (x)[H,x l ](j)j(x) = i(E n - Ej) J dx (j)\ l (x)x l (j)j{x) . (40) 
The correction due to variation of the bound state eigenvalue is thus 

5 E £ L = -5E n + -P / dz5 E U(z) , (41) 

7T 71 JO 

where 

poo poo 2 

5e14(z) = —5E n / dx 2 x\ / dx\ x\ V" V" 

JO JO „ 



K i,j=l 

d 

h(x 2 )GV(x 2 ,x ll z)fo(x l )—A l *{x 2 ,x l ,z) . (42) 

From Eq. ([I]) one can also easily obtain the correction due to the variation of the bound 
state wave function: 

2a T? f En j_ f j_ fj„ j.U-\~.i^,- ~ ^„.m- / x/r « w™sin[(£ n -z)a; 2 i] 



5 <A = —^ F J " dz J dx * J dx ± <t>i(x 2 )a l G(x 2 , x x , z)a m 5(j) n (x 1 )(5i m V 2 ■ Vi - V 2 V?" 

since the dependence on the bound-state wave function is explicit. 

D. Low-Energy Green's function correction 

The correction due to the variation of the Green's function in Eq. flT7|) is given by 



(E n - z) 2 x 21 
(43) 



6 G S L = -Rel dzd G U(z) (44) 



7T JC 4 

where first-order change in 14 (z) in Eq. fll8|), due to variation of Q is 



/•oo poo 

5qU{z) = — / c?X2 ^2 / dx\x\ 
Jo Jo 



2 

E fn,i(x 2 )5G l i(x 2 ,x 1 ,z)f nJ (x 1 )A l 3(x 2 ,x 1 ,z) , 

(45) 

and the integration contour C + is defined subsequently. The Green's function Q l J(x 2 , x±, z) 
in Eq. ([Tj|) has poles along the real axis in the range of integration over z, and in the 
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perturbation expansion of G^(x 2 ,xi, z), in powers of a perturbing potential 5V and the 
unperturbed Green's function G]j>(x 2 ,Xi,z), higher order poles are introduced as noted in 
Sec. [TI|. Here, our method of isolating those poles and evaluating their contribution is 
described. In terms of the spectral resolution of the unperturbed radial Green's function 

fm,i\ x V fm,j( x l) 



(46) 



Eq. flTED reads 



/■oo 2 

SG^(x 2 ,x 1} z) = - dx 3 x 2 3 J2 E 

k=l rn 2 ,m 1 

fm2,i\. x 2)fm2,k\ x 3) fry ( x ^ fmi,k\ x 3) fmi ,j\ x l) 



E m2 z 



(47) 



which explicitly shows the second and first-order poles. In (|7), only states rri2,mi with 
spin-angular momentum quantum n contribute. The principal parts of 5G^(x 2 , Xi, z) are 



identified by expanding the functions in Eq. (f4"5|) in Laurent series about each pole. For 
z ~ E m , 

frn,i( x 2)fm,j( x l) 



G%(x 2 ,xi,z) 



E - z 

J -'m " 



G^(x 2 ,x l ,E mJ 



+0{z-E m ), (48) 

where G^i (x 2 ,x\, E m ) are the radial components of the reduced Green's function given in 
Eq. © 



G^' ij (x 2 ,x 1 ,E n 



lim 

Z—*E„, 



E 

i 



G^(x 2 ,x 1 ,z) 

frriji^X^fmJ 



Em z 

fl,i( X 2)fl,jM 
El — Em. 



(49) 



Hence 

5G^(x 2 ,x 1 ,z) 



r 1 

dx 3 X% f m A X 2)frn,k{x^)5V(x Z )f mM (x^)f m j(Xi)— — 

k=1 L {& m - z) 



+G^ lk (x 2l x 3 , E m )5V(x 3 )f m:k (x 3 )f m j(x 
+fm,i(x 2 )f mjk ( y x 3 )5V(x 3 )G^ kj ( y x 3 , Xi, E n 



1 



E„ 



Em, Z 



0{1) 
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or, in view of Eqs. @ and 



SG%(X 2 ,XI,Z) = -f m ,i{x 2 )f m j{x X )- 



5E r , 



+ 



8fm,i(x 2 ) fm,j( x l) + fm,i{x 2 ) ^fm,j{X\j 



1 



Ern % 



(E m - zf 

The same result may be obtained by expanding the pole contribution to the full Green's 
function 



+ 0(1) 



G K (x 2 ,Xi,z) 
in powers of SV(x), with 



&m. Z 



&m E m + 5E m + • ■ ■ 



and retaining only the first-order correction 

~fm,i( x 2)fm,j( x l) 

+ 



(E m - zf 

5fm,ii X 2) fm,j( x l) + fm,i( X 2) $fm,j( x l) 



E - z 



In addition to the expansion of the Green's function correction, we have 

A%(x 2 , x u z)= A i ^(x 2 , x u E m ) 

+{z-E m )B i i{x 2 ,x 1 ,E m ) 



where 



+Ol{z-E n 



B i K j {x 2 ,x 1 ,E m ) = — [A%{x2,xx,z) 



z=E„ 



The complete expansion for z ~ E m is 



dx 3 x 2 3 G^(x 2 ,x 3 ,z)8V(x 3 )G k K 3 (x 3 ,x 1 ,z) A^(x 2 ,xx,z 

1 



k=l 

fm,i{ x %) frn,j( x l) 5E„ 



A^(x 2 ,x u E n 



B^(x 2 ,x u E v 



( E m z)^ E m z 

Sfm,i( X 2) fm,j( x l) + fm,i{ X 2) S fm,j i x l) A l ^{x i ,X 1 ,E, 



(52) 



(53) 



(54) 



(55) 



Ern, % 



0(1) . (56) 
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Our strategy for dealing with poles in the low-energy part is to calculate the line integral 
over z of the difference between the complete integrand and the pole terms, and add the pole 
terms integrated analytically. In the unperturbed self-energy calculation, the singularities 
along the real axis in the interval (0, 1) are poles, and the appropriate prescription for 
integration over z yields the principal value integral in Eq. flT7|). In the present context, 
there are double poles as well, so it is necessary to reexamine the original derivation to 
obtain the correct prescription. It follows from the discussion in Ref. || that the integration 
over z can be written as 



Re / dzU(z) 
Jc+ 



(57) 



where C+ is a contour that extends from z = to z = E n above the real axis in the complex z 
plane. Here we use a method based on an analytic evaluation of the pole terms, as described 



in Ref. | 14fl . With the notation IA V for the pole terms of U, we have 



Re f dzU(z) = Re f dz [U(z) -U p (z) +U p (z) 



rE n 

/ dz 
Jo 



U{z)-UJz) 



1 



—Re ^ dz'U p {z') 



(58) 



The pole at z = E n , the endpoint of the integral over z, does not cause any problem, as 
follows from the discussion in Sec. |ill (J| [see Eqs. ([57]) to (pUf)]. The relevant integrals for 
the analytic evaluation of the pole terms are 

1 E n 



Re / dz 



c+ (E m - zf 



E m (E n E m j 



(59) 



where 



Re / dz- 



c+ E m z 



In 



E n E m 



(60) 



< E m < E n 



(61) 



Applying these results to Eq. (0), we have 
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a f 
5g£l = — Re / dzbo,lA (z) 
7T Jc + 

/ \8 G U{z)- £ 
Jo ™ 



a r En 

71 



m 

E m <E n 



II .III ^ ' II .III 



(E m - zf (E m - z) 



E m <E n 



F?W E / E 



E m {.E n E ri 



E - E 

J -'n J -'m 



(62) 



where 



i?l 2 L = SE m / rfx - 2 
/o 



oo 



f°° , - 

„.,„ - <<< , / ..-2 2^ / dXiX 2 £ / ni ^ 2 )/ m ,i(x 2 )/ m j(^l)/nj(^l)^«(^2,^l,-Em) (63) 



and 



/"°° /"°° 2 

R nm = -S E ™ / ^2^2 / ^1 ^1 H fn,i{x 2 ) f m ,i (x 2 ) fm ,j (%1 ) fn j(Xl ) -B* 3 (x 2 , X\ , -E m ) 
•/() •/() 



2 x 2 i dx x x\ J2 fn,t(x 2 )[5 f m:i (x 2 ) f m ,j(xx) + f mii {x 2 )8f m j{xi)\f n j{xi)A^{x 2 ,x l ,E n 
i,j=l 



dxr> x 2 



E. Numerical evaluation of the first-order correction to the Green's function for the 

low-energy part 



Numerical evaluation of the Coulomb Green's functions in this paper is based on the 
explicit formulas given, for example, in Eqs. (A. 16) and (A. 17) of Ref. 0, together with 
the numerical algorithms described in Ref. ||15|| . The first-order correction to the Green's 
function in Eq. (|16|), for the range of arguments relevant to the low-energy part, is evaluated 
by numerical integration over x%, where the interval of integration is divided into four subin- 



tervals: (0,yi), (yi,y 2 ), (2/2,2/3), and (y 3 , 00). Defining p = 2y/l - z 2 , y = 3, and assuming 
x 2 < xi, we choose y x = mm(y ,px 2 ) } y 2 = m.in(ma.x(y ,px 2 ),pxi) } and y 3 = max(y Q ,pxi). 
In the first interval we make the substitution £3 = £ 3 , and integrate over £ by Gauss-Legendre 
quadrature with 12 to 26 integration points. In the second and third interval we also do 
Gauss-Legendre integration, with 17 to 33 and 11 to 21 points, respectively. We use 6 to 
18 point Gauss-Laguerre integration for the remaining interval. The integrations over the 
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second and fourth intervals are the least accurate. For a 2s electron at Z = 20 the integral 
over (7/1,7/2) has an error of a few parts in 10 6 in the worst case. 



F. Reduced Green's Function 



In the preliminary version of this calculation ||10|| , a purely numerical method of evalu- 
ating the reduced Green's function was employed. However, while that method is adequate 
at high Z, it gives unsatisfactory results for the 2s state when Z < 20, so a new method 
was developed that yields better precision. As a check of the coding of the later method, 
the results of the two methods were compared and are in agreement within a relative differ- 
ence of 10 -6 over a wide range of the variables. Both methods are briefly described in the 
following subsections. 



1. Evaluation of the Reduced Green's Function by numerical pole removal 
Eq. (|15D can be written as 

G«(x 2 ,x 1 ,E n + e) = - / "'* (X2) e / "' j(Xl) 

+Gy j (x 2 ,x l ,E n ) + 0(e) . (65) 

As an immediate consequence of this relation, we have 

^[G%(x 2 , xi, E n + e) + G*(x 2 , x 1 ,E n -e)] 

= G^(x 2 ,x u E n ) + 0(e 2 ) , (66) 

so the reduced Green's function can be easily obtained from the full Green's function by 
symmetric interpolation of the energy variable. We form a linear combination of two such 
interpolations in order to obtain a result with an error of order e 4 rather than e 2 . In 
particular, we have 

|[G^'(x 2 , x h E n + e)+ GH(x 2 , xi, E n - e)] 
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[G*(x 2 , x u E n + e') + G^(x 2 , x u E n - e')] 



G*"(z 2 , x l5 E n ) + C 2 [ge 2 + (1 - q) e' 2 ' 
+C 4 fge 4 + (1 - q) e ,4 l + 0(e 6 ) + 0(e' 6 



(67) 



where the choice 



g > 1 



(68) 



provides 



qe 2 + (l-q) e' 2 = 



(69) 



with the parameter q free to vary in the range (1, oo). The choice q — | gives the Lagrange 
interpolation formula with equally spaced evaluation points: e' = 2e. With the choice in 
Eq. (|68|), the fourth-order term is proportional to 



ge 4 + (l-g)e 



9 



(70) 



The coefficient of e 4 is 4 for equally spaced points and approaches 1 as q —>■ oo. As a 
compromise between a large coefficient for q near 1 and the minimum coefficient as q — > oo, 
with a correspondingly larger roundoff error, we employ the value q = 10. The interpolation 
interval is taken to be e = 0.02 (E n+ i — E n ) where E n and E n+ i are the Dirac eigenvalues 
for principal quantum number n and n + 1 for the same k. This interval avoids overlap with 
the nearest pole of the Green's function. 



2. Direct evaluation of the Reduced Green's Function 



The method developed for the present work is similar to that of Hylton ||T3|,|16j, but 



differs in the details of its implementation. From Eq. (pop, it is evident that we obtain the 
reduced Green's function by expanding the various components in the explicit expression for 
the radial Green's function in powers of e = z — E n and keeping only the final combinations 
of terms that are of order 1. To implement this, at each step in the numerical evaluation 

17 



of the reduced Green's function we calculate only the coefficients of the leading two terms 
in the power series in e and discard the higher-order terms. In certain cases, it is necessary 
to begin with three terms in the expansions, because the leading term either vanishes or 
cancels an equal leading term in forming a difference. As suggested by the form of the 
following equations, many of the coefficients follow from combinations of coefficients that 
appear earlier in the calculation. 

The code for the numerical calculation was written by modifying the existing code for the 
radial Green's function described in Refs. and [|l7j], so only a few details that illustrate 



the approach are given here. We define the radial quantum number n r = n—\n\. Expansions 
are needed for z: 



z = z + z X 6 + z 2 e 2 H 

z = E n ; z x = 1; z 2 = , (71) 



for c = y/1 — z 2 : 



c = c + cie + c 2 e 2 + ■ 



7-E n n r + A 1 

Co = ^Ta ; Cl = -— ] C2 = "M' (72) 



where 7 = Za and A = ^k' 2 — 7 2 , for w = j/c: 



w = w + Wit + w 2 e 2 H 

7 7Ci l{c\ - c c 2 ) 

w = — ; wi = r ; w 2 = 3 , (73) 

c c c 



and for v = •yz/c: 



v = u + v x e + z/ 2 e 2 H 

v = z w ; v x = z w 1 + z^o; 

v 2 = z w 2 + ZiWx . (74) 
Expressions that appear in the definitions of the radial Green's functions include A — v. 
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A — v = —n r — v\€ — i>2e -\ 

| — v\e — u 2 e 2 + • • • if n r =0 
—n T — v\E + • • • if n T ^ 



and r(X — v): 



r(A - 1/) = r^V 1 + r w + 



Z/i z/f 
Ji("r) ]_ p(n t -i) 

r (n r ) = _l r K-D + ^ r K-D ( 76 ) 

where 7e = 0.577 ... is Euler's constant. The recursion relations used to calculate the power 
series for the Whittaker functions are treated in a similar manner. For example, for the 
power series evaluation of M v+ (i/2),\(2cx) by means of the recursion relations in Eq. (D.2) 
of [|15|], we write 

T{m) = T (m) + Ti(m)e + • • • 

T (0) = 1; T 1 (0)=0 

m / \ (m — n r )2cnx ^ , . 

T (m + 1 = t — : — --— -T m 

0V ; (m + 2A + l)(m+ 1) 0V ; 

^ , \ (m — n T )2c x ^ , . 

TAm + 1 = i— : -Tim 

V ' (m + 2A + l)(m + 1) u ; 

i [(m-n r )c 1 -, lCo ]2x To(m) (??) 



(m + 2A + l)(m + 1) 

The termination of the power series for the leading term in Eq. (|77|) at m = n r corresponds 
to fact that the leading term in Eq. ( |65l) is proportional to the bound-state wave function. 

The calculation of the reduced Green's function in this work is based on the application 
of two-term expansions in e, as described above, to the numerical evaluation of the complete 



Green's function (see Appendix D of Ref. The numerical value of the reduced Green's 



function is just the collection of terms with combined order 1 in e. 
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G. Numerical evaluation of the first-order correction to the wave function 



The first-order correction to the wave function, given by Eq. (12), is evaluated with the 
aid of the reduced Green's function as described in Sec. ffTT]. The numerical integration in 



Eq. (12) is divided into 3 segments: 

Pn 



f 1 I \ 2 

5 / duxiu A (xtu 5 /p n j ^2G^ v (x 1 x 1 u 5 /p n ,E n )6V(x 1 u 5 /p n )f n>j (x 1 u 5 /p r , 

J ° 3=1 

+ r dy (y/p n f G™(x, y/p n , E n )5V(y /p n )f nJ (y/p n ) 



3=1 

+ / dy(y/p n fJ2 G ? ij (x,y/Pn,E n )8V(y/p n )f n!j (y/p n ) , (78) 



where p n = 2y 1 — E% is the coefficient of the argument in the exponent that governs the 
behavior of the integrand for large values of the argument, x\ = min(xp, 3), and x<i = 
max (xp, 3). The first and second integrals are evaluated by means of 30 point Gauss- 
Legendre quadrature and the third is evaluated with 15 point Gauss-Laguerre quadrature. 



The accuracy of the perturbed wave function calculated according to Eq. (|78|) has been tested 
in the Coulomb case by comparison to the result obtained by numerical differentiation of 
the Coulomb wave function. 

For values of the argument of the first-order correction to the wave function near the 
origin, we found that greater numerical accuracy and speed could be obtained with a nu- 
merical evaluation based on the expansion in powers of x. This expansion is described in 
Appendix [B]. 

IV. HIGH ENERGY PART 

The high-energy part, given by the integral in Eq. (|^), must be regularized, since it is 
formally infinite. We employ the Pauli-Villars regularization scheme, following the method of 
Refs. [1(],12] to isolate and remove the divergent contributions. In Ref. [12] we demonstrated 



that suitable numerical convergence can be achieved through the use of a term-by-term 
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subtraction method. This particular method has the advantage that it does not require a 
mix of coordinate-space and momentum-space calculations as do earlier methods, but works 
entirely within coordinate space. In this method, the high-energy part £h, given by Eq. (0), 
is separated into two parts: £ha and £hb- The divergences are all contained in £ha and can 
be calculated completely analytically, while £hb is finite and is treated numerically. In this 
section we describe the method used to evaluate while the method to compute 5 Era 

is discussed in Sec. [V|; the total is 



SEn = 5E HA + SE HB . 



(79) 



The high-energy remainder with term-by-term subtraction, from Eq. (32) in f[2" |, is 
written as 



£hb — £h — £ha 

a 



oo poo 

2 J™ J2 



, dz / dx 2 x 2 / dxi x 1 
2tti Jc h Jo Jo 



< ^ K (x 2 ,x 1 ,z)-lC^\x 2 ,x 1 ,z)-lC^ 1 \x 2 ,x 1 ,z)-]C^ \x 2 ,x 1 , 
|#c|=i 



K,§' 2 \x 2 , Xl ,z) \ , (I 



with 



K. K (x 2 ,xi,z) = [^ 7 n,i( x 2)SK( x 2,xi,z)J : ' n j(xi)A K (x 2 ,xi) - T n ^{x 2 )Q^ {x 2: x x , z)T nJ {x 1 )A^{x 2: x x )\ , (8 
»',i=i 



JC^°\x 2 ,x u z) = A K {E 1 K \x 2 ,x 1 ,z) [J* tl (x 2 ) -3Fl 2 (x 2 )} +E 2 l(x 2 ,x l ,z) [j*j(x 2 ) - 3J* tl (x 2 )]} , (82 



= {f?(x 2 ,x 1 ,z) 
+F 22 K (x 2 ,x 1 ,z) 



F n ,i(x 2 ) (b k -^- + C k 



F n , 2 {x 2 ) B K 



_d_ 

dx 2 
d 



C K 



+2 Fii(x 2} Xl , z)T n>1 {x 2 ) \ B K — + C K 

d 



+2F% 1 (x 2 ,x 1 ,z)F n , 2 (x 2 ) ^B K ^ + cJ-^j ^lO**)} 



1 






x 2 


1 








1 






^2 


1 


+ fv n 



-^,1(^2) - ^,2(^2) \ $BK-jr — 



' dx 2 

d 



Tn, 2 {x 2 ) - T n ,\{x 2 ) 3B K - C K 



x 2 , 

1 + K n 



dx 2 



x 2 



Fn,2(x 2 ) 
I Fn,l( X 2) 



Fn, 2 (x 2 ) 
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and 



d 



JC^(x 2 ,x 1 ,z) = -V(x 2 )Aj F 1 K \x 2 ,x 1 ,z) J* tl (x 2 ) - 3^ n 2 2 (x 2 ) 



dz 



(84) 



In Eqs. ( PH ) to (Q) A^, A K , _B K and C K are integrals over coordinate directions denned in 
Refs. P, [l2|,[l4],|l8l , T n> i{x) are radial components of the bound-state wave function as before, 
and K n is the spin-angular momentum quantum number of the bound state n. The ex- 
pressions for the free Green's function radial components F^(x 2 , X\, z) and their derivatives 
can be found in |12[| and those of the Coulomb Green's function can be found in [Ell. The 



methods we used for summation over angular momentum k and for numerical integrations 
are identical to those described in [[TJ] and will not be repeated here. We also found that 
the convergence of the numerical integration was much better than in the case of the un- 
perturbed self-energy. We thus did not use the extra subtraction term K^' 2 \x 2 , 21, z) that 



was necessary in [12| to obtain good convergence at low Z. The high-energy remainder 
for the self-energy screening is obtained from Eq. fl8"0D as described in Sec. |TJ as the sum 
<^£hb + ^g^hb + <5e£hb- In the three following subsections we derive the expressions that 
are used to obtain ^^hb, <5g£hb, and <5e£hb from Eqs. (|80|) to 



A. Wave function correction 



To obtain the expression for the high-energy remainder for the wave function correction, 
we need the functional derivatives of Eqs. (|8~1~D to ( |84] ) with respect to the radial wave 
functions T n ^{x). For the full expression we have 
2 

K^ K (x 2 ,Xi,z) = { [fn,i( x 2)G^(x 2 ,x 1 ,z) 5f nd (xi) + Sf nji (x 2 )G^(x 2 ,x 1 ,z)f nd (xi) A K (x 2 ,Xi) 



fn,i(x 2 )G^(x 2} x 1} z)5f nJ (x 1 ) + 8f n>% (x 2 )G^ '(x 2 , x x , z)/ nJ (:ci) A^(x 2 ,Xi)} , 



and for the subtraction terms, we obtain 
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K f,'° \ x 2,Xi,z) = 2A k {F^ 1 (x 2 ,X 1 ,z) [/„,i(x 2 ) Sfn :1 (x 2 ) -SfntMSfnjfa)] 

+F* 2 K (x 2 , xi, z) [f n ,2{x 2 ) 5f Ht2 (x 2 ) - 3f n>l (x 2 ) 8 f n>l (x 2 )} } , (86) 



Kfi\x 2 ,x X: z) 



F K (x 2 ,x x ,z) 



fn,l{x 2 ) \ B R -^— + C K 1 + Kn ) 8f nA (x 2 ) - f n ,2{x 2 ) \3B K -^- C K - —J 8 f n , 2 {x 2 ) 



"dxo 



d 



+Sfn,lM [B K — + C K 



x 2 

1 + K n 



dxo 



d 



x 2 

1 K~ 



fn,l( X z) - S fn,2( X l) 3B K - C K 

X 2 I V " x 2 X 2 



fn,2{x 2 ) 



+Ff 2 K (x 2 ,x 1 ,z) 



fn,2(x 2 ) ( B K ^— + C K 1 8f n ,2{x 2 ) - fn,l( X 2) ( 35 «T^ ~ C K 1 + — ) 5/ n ,l(x 2 ) 



(9x9 



<9 



+8f n)2 {x 2 ) \B R — + C K 



x 2 

1 — K r 



dx-2 







x 2 



fn,2(x 2 ) - 8f nil (x 2 ) \ 3B K — C K 



x 2 

1 + K r 



X 2 



fn,l{x 2 ) 



+2F! 2 K (x 2 ,x u z] 
+2F* 1 (x 2 ,x h z) 



fn,l(x 2 ) (B K ^- + Cj^^\ 8f n , 2 (x 2 ) + 8f n ^X 2 ) (b k ^- + Cj^^\ fn,2(x 2 ) 







1 + K. 



fn,2{x 2 )\B K —+C K 

\ ox 2 x 2 



8fn,l{x 2 ) + 8f n<2 {x 2 ) B K 



_d_ 

dxo 



C K 



1 + Kn 

X 2 , 



fn,l{x 2 ) 



and 



K $,k\ x 2i x 1i z ) = -^y{ x 2)A K \^—Fl 1 (x 2 ,Xi,z) [f n ,l(x 2 ) 8 f nA (x 2 ) - 3f n>2 (x 2 ) 8fn i2 {x 2 )\ 

+^F 22 K (x 2 , x x , z) [fn,2(x 2 )8f ni2 (x 2 ) -3f ni i(x 2 )8f nA (x 2 )] 



In terms of the expressions in Eqs. fl55|) to (p8|), the first-order wave function correction to 
£hb is 



q. p roo poo 

8(h£iiB = dz dx 2 x\ \ dx\ x\ 
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x K <i>A x ^ x ^ z ) - K l,n ( x 2, x i, z ) ~ K^' K '(x^x^z) - Kly[x % ,x u 

\k\=X 



39) 



In order to evaluate the expression in fl3?D, we need the derivative of the bound-state 
Dirac wave function and the derivative of its first-order correction in the potential 8V . The 
differential equations for the large and small components of the unperturbed wave function 
are (see, e.g., |TD|, Appendix A) 

^fn,l( X ) = -^^fnAx) + [1 + E n - V{X)\ f n , 2 (x) 
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= [1 - E n + V(x)) f n>1 (x) - L_*»/ ( x ) , 

ax x 



(90) 



which yield the wave function derivatives from the analytic expressions for the wave function. 
We obtain analogous expressions for the perturbation of the wave-function components in 
the potential V(x) + 5V(x). Retaining only first-order terms in SV(x), we obtain 

^/n,i(x) = -^-^Sfn,i(x) + [1 + E n - V(x)} 6f n>2 (x) + [6E n - 6V(x)]f n , 2 (x) 
^k,2{x) = [1 - E n + V(x)} 5f n> i(x) - ^—^8f n , 2 {x) + [5V(x) - 5E n ]f n> i(x) . (91) 



B. Green's function correction 



We do the corresponding calculation for the high-energy term to account for variation of 
the Coulomb Green's function under a change of the potential. Since the Coulomb Green's 
function has no poles on the high energy integration contour (which lies on the imaginary 



axis), we may directly apply Eq. ( I6[) to obtain 



2 /.oo r 

K G:K (x 2l x u z) = - / dx 3 x\5V{x 3 ) f n , i (x 2 )G lk (x 2 ,x 3j z)G k . : >(x 3j x 1 ,z)f n j(x 1 )A K (x 2 ,x 1 ] 



i,j,k=l ' 



-f n j(x 2 )Gf(x 2 ,x 3 ,z)G k K j (x 3 ,x 1 ,z)f nJ (x 1 )A i ^x 2 ,x 1 ) 



(92) 



Only f^ 1 ' ) contributes to the subtraction term, and we thus obtain from Eq. 



K i G^\x 2 ,x l ,z) 



-8V{x 2 )A K \ ±Fl 1 (x 2 ,x 1 ,z) [fl^) - 3/ 2 2 (* 2 )] + j-F^x^z) [/ n 2 2 (x 2 ) - VlM) 



The Green's function correction to £hb is thus 



5 a E 



G<--HB 



— / dz dx 2 x\\ dxi xj V] K Gk (x 2 ,x 1 ,z) - > {x 2 ,x 1 , z 
2ni Jc H Jo Jo ,^z 1 L 



(94) 
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C. Numerical evaluation of the first-order correction to the Green's function for the 

high-energy part 



The evaluation of the the Green's function correction to the high-energy part is the 
most difficult and time consuming part of the present calculation because of the additional 
non-trivial integration over x 3 and the necessity of evaluating integrals over a product of 
two Coulomb Green's functions in Eq. fl92|) . The calculation is facilitated somewhat by 
subtracting and adding an additional term that closely approximates the term in Eq. fl92|). 
We write 



dx 3 xl 5V(x 3 ) G ik (x 2 , x 3 , z) G kj (x 3 , x u z) 



dx 3 xl 



5V(x 3 ) - l -5V{ Xl ) - l -5V{x 2 ) 



Gf{x 2 ,x 3 ,z) G^{x 3 ,x u z) 



where (see Eq. (35) from Ref. fllPl ) 



+-[5V( Xl ) + 5V(x 2 )\ / dx 3 x 2 3 G ik (x 2l x 37 z)G k J(x 3 ,x u z) 
2 Jo 



V / dx 3 x 2 3 G l K {x 2 ,x 3 ,z)G k3 {x 3 ,Xi,z) = —G %3 (x 2 ,x 1 ,z) . 



We calculate separately the two terms in the sum 5g^hb = ^g^hb + ^g^hb' where 



'(2) 



m en r f°° f°° 00 ( \ 2 

^g^hb = 7T-- \ dz dx 2 x\\ dx\x\ Y] \ - - V] [5V(xi) + 
2tti Jc n Jo Jo |k|=1 I 2 iJ=1 



SV(x 2 )} 



x 



(95) 



(96) 



d i) ... , 

fnAX2)g-G^(x 2 ,x 1 ,z)f n j(xi)A K (x 2 ,x 1 ) - f n - l (x 2 )—G l3 (x 2 ,x 1 ,z)f nJ (x 1 )A t3 (x 2 ,x 1 ) 



(97) 



and 



x c(2) 
^g^-hb 



a 



OO /"OO 

dx 2 xl 



poo rc 

I dxi xl 
Jo 



dx 3 xl 



SV{x 3 ) + Uv{ Xl ) + Uv{x 2 ) 



dz 

2m Jc u Jo 

oo 2 

x [fn,i(x 2 )G lk (x 2 ,x 3 ,z)G k3 (x 3 ,x 1 ,z)f nJ (x 1 )A K (x 2 ,x 1 

[re|=l i,j,k=l 



-fn,^x 2 ) G ik (x 2} x 3 , z) G kj (x 3 , xi, z)f nJ (x 1 )A i i(x 2 , x x )] ■ 



(98) 



The derivative of the Green's function in Eq. fl9"T| ) is evaluated numerically with a two-point 
formula: 
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D. Energy correction 



We evaluate 5e£hb, given by 

<5e£hb 



SE 







n dS r 



■s 



HB 



(100) 



by numerical differentiation of £hb- We employ an effective infinitesimal displacement of 
the energy variable obtained by making the replacement of E n (Za) by E n \{Z ± e)a], with 
e = 0.002. The symmetric derivative formula with this displacement has an uncertainty of 



a few parts in 10 . One expects that numerical integration errors, as discussed in [12 



arc 



slowly varying functions of Z and largely cancel in forming the derivative, leading to an 
accurate result for the derivative. Evidently, however, the error in terminating the sum over 
k does not cancel as strongly, which may lead to significant uncertainties for some Z and n£. 
We have thus found it preferable to differentiate directly the running term in the sum over 
angular momenta in Eq. (|80|). To study this effect, we compare numerically -^-Sa{t, y, u) 
evaluated by differentiation of the sum 







dE 7 



■SA(r,y,u) 







E 



dE n 

-Kff ( r V, V, iu) - K§£> (ry, y, iu) 
-K^f\ry,y,iu) 



(101) 



to a sum of derivatives (term-by-term differentiation) 

K^ K {ry,y,iu) 



dE. 



d 00 d 

S A {r,y,u) = J2 



K \=l 8E n 



- K ff (?!/, y, iu) - K { ^ ] {ry, y, iu) 

-K^firy^y,^) 



(102) 



where u = (1/2) (1/t - t) and are defined in Eqs. 



to 



. The summation is 



terminated when the remainder estimate as described in Ref. [12| is smaller than a predefined 
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cutoff value. Although the two methods of calculation converge to the same value for a very 
small cutoff (10~ 10 ), the results with a larger cutoff can disagree by a few parts in 10~ 4 . 
Since very small cutoff values would lead to other problems when r is close to 1, we evaluate 
the energy derivative by the term-by-term method. Other illustrations of this problem are 
discussed in Sec. |V1 A| . 



V. ANALYTIC TERMS 

In this section we evaluate analytically the terms that are subtracted in the numerical 
calculation described in Sec. |IV| for each of the three contributions to the screened self-energy. 



In Refs. flQj]12]l , with this application in mind, we were careful not to obscure the origins 
of V and E n through the use of the Dirac equation. The individual terms corresponding to 
the three diagrams of Fig. are thus obtained by direct differentiation. 



A. Wave- function correction 



The analytic portion of the wave-function correction follows from Refs. ||10yi2|| . In each 
of the terms E# , we calculate the variation with respect to a change in the wave function 
based on the explicit dependence on the wave function, and obtain 

2a 



(0,0) 



7T 



: (</>n|/Wn> 



hi(^) - 1 + l —J^L In (l + El) + 0{A- 1 ) 

-*—Jrn 



(103) 



(0,1) 



and 



— (<p n \at-p \5(p n ) 

7T 



6 - 3E? + 7Et 1 - Ei 



2AEl{l + El) AE^ 



In (l + El) + OiA- 1 ) 



(104) 



Srh£ 



(1,0) 



2a 



71 



(<i>n\V\6<i> n ) 



6 -El 3 + Ei 



AEi 



- In (l + El) 



(<j> n \pv\6<f> n ) 



2 2 , 
E--El l < l + E ^ 



+ 0{A 



(105) 
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In analogy with [|T(],|T^], we define 



(^£ha = lim 



+5^ 0) -2 5m(A) (0 n |/3|<50 n > 



(106) 



where the last term is the renormalization term and 



ol 

5m(A) = — 

7T 



~lnM 2 ) + - 
4 K ' 8 



(107) 



If we combine the coefficients of ln(yl 2 ) in Eq. (|106|) , we obtain 



a a 
— (<p n \(3 + a ■ p + V \6<f> n ) = —E n (<j) n \5<j) n ) = 
2tt Lix 



(108) 



from the differential equation for <ft n and the fact that the first-order correction to the 
wave function is orthogonal to the unperturbed wave function. Hence, the wave function 
correction is separately finite. 

B. Energy correction 

This correction is obtained by differentiating all terms with respect to S n . The three 
terms are 



x c(o,o) 



a 



SE n 



71 



2 I- El 



Kl + El El 



In (l + E. 



6-E 2 3 + E, 



+ *ln l + £ 2 )+0(A- 1 



(109) 



s c(0,l) 



and 



71 



5E n (<f> n \a ■ p\(j) n ) 



6 + 9E 2 - 8E* - 3Ef, 1 



6^(1 + ^2)2 El 



log(l + El) + OiA- 1 ) 



(110) 



6 E £g' 0) = -SE n \ {<t> n \V\<t> n ) 

71 



6 + 3EI + E* 3 , M 

logl + £ 2 



+ {<t>nW\<t> n ) 



6 + 2El 



m(\ I E; E\\ 
6 



E 2 (l + E 2 ) E 4 



log(l + # 



+ e>(yT 1 )j . (ill) 
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The total is 



#e£ha = <5e£h°'°' ) + ^E^n' 1 ' + ^E^H 



(0,1) 



,(1,0) 



(112) 



which contains a divergent term given by 

a 5E n 



e<-ha 



An 



ln(A 2 ) + 0(1) . 



113) 



C. Green's function correction 

This term arises entirely from , which is linear in the potential in the Green's 
function. Taking into account the fact that {<j) n \5V\<j) n ) = SE n , we obtain 



5 G £ RA = 5 G £S' 0) = -{5E P 

7T 



1 in(^) + ^-^m(i + ^) 



8EZ 



-(<j> n \psv\<i> n ) 

This expression also contains a divergent contribution 

<5g£ha 



+ 0(A 



-i^ 



a6En HA') + o(i) 



47T 



which cancels the corresponding term in Eq. ( |113| ). The total expression 



(114) 



(115) 



(116) 



is thus finite as expected. 



VI. RESULTS AND DISCUSSION 
A. Coulomb tests 

In order to check the accuracy of the numerical calculation, the equations derived in 
the previous sections, and the parallelized code, we compare the calculations obtained by 
numerical differentiation of the one-electron self-energy function F(Za) to the results of 
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the method presented in this paper using 5V(r) = —ea/r as a perturbing potential. This 
method, proposed in Ref. is very efficient, as each individual contribution to the screened 
self-energy can be checked independently. The overall agreement between the results of these 
two methods of calculation is good, although differences between some contributions can be 
several times the combined uncertainties based only on the apparent convergence of the 
numerical integration. These additional errors come from the numerical problems described 
in Sec. |IV D| . An illustrative example is the wave function correction. One can compare the 



term-by term derivative with respect to Z in the wave function to the derivative obtained 
from two converged sums, as is done in Eqs. ( |101| ) and (|102|) . There is also an independent 
calculation based on the first-order correction to the wave function from Eq. @. The 
evolution of this sum for smaller and smaller values of the cutoff error are displayed in 
Fig. U One can see that although the two calculations converge to the same limit when the 
cutoff is as small as 10 -15 , the results follow very different paths. Only the term-by-term 
differentiation method follows the result obtained from evaluation of the first-order correction 
to the wave function independently of the cutoff. This constitutes a very demanding test of 
our numerical evaluation of the first-order correction to the wave function. The difference 
between the two calculations is never smaller than 8 x 10 8 , which is the error from numerical 
uncertainties of the full Green's function and the error in the numerical derivative. 

To improve the numerical precision, we have employed a subdivision of the integration 



over r into regions with < r < 0.4 and with 0.4 < r < 1, as described in Ref. ||12|| . This 
division provides an accurate evaluation that does not require functional evaluations with 
values of r too close to 1. In this way we where able to obtain an accurate comparison of 
all contributions in the high-energy part. 

A few such problems remain in the test calculations in the low-energy part at low Z. We 
did not attempt to improve the accuracy, because we have enough accurate cases to check 
the code and numerical procedures. From the tests we have performed, it is clear that the 
calculation based on numerical differentiation is the less accurate. However since numerical 
differentiation is used in the final result to obtain the reducible correction we have increased 
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the total uncertainty of both the pure Coulomb test and spherically-averaged potential of 
the next section accordingly. 



B. Results with spherically-averaged one-electron potential 

The calculations of interest for physical applications are based on realistic potentials 
obtained from the spherically-averaged potential of Eq. flC^ ) of Appendix [C[ All results 
presented here are given in terms of the scaled function defined by 

£n£j,n'e'j' — — ^ 3 ' Fnej,n'e'j>(Zo>) mc 2 (117) 

All 16 possible total scaled functions FiZct) for the self-energy screening of ni electron 
by a mil' electron, 1 < n, n' < 2, < £, £' < 1 are given in Table and in Figs. |5] to ^. It can 
be seen that the uncertainty at low Z can be as high as 30% for the screening of 2p electrons 
at Z — 10, or as low as 10 -6 . These functions can be used to evaluate the self-energy 
screening correction to any atom with two to ten electrons, in a shell n < 2. As an example 
we treat the case of lithiumlike uranium. With the results presented here we can compute 
the self-energy screening correction for all three states ls 2 2s, ls 2 2pi/ 2 and ls 2 2p 3 / 2 . A first 
approximation is obtained for transition energies by neglecting the core relaxation. The 
self-energy screening correction to 2pj — > 2s transition energy is evaluated using Eq. (|H7|) 
as ASj = 2^2pj,is — 2^2s,is, where the factor of two accounts for the fact that there are two 
Is electrons screening n = 2 electron. A better approximation, which takes into account the 
relaxation of the core electrons, and provides a value for the total binding energy is given 
by Si s 2 2s = 2(£ 2S) i s + £i St 2s + £u,u) and £is 2 2 Pj = 2(£ 2 p 3 ,is + £u,2 Pj + £i s ,i s )- The results of 



these calculations are presented in Table |V| together with all other calculations known to 
date. 

It should be noted that the present method is equivalent to the Coulomb approximation 
in certain cases. If one considers only the Coulomb contribution in the interaction between 
the two electrons in Fig. 0, the present method provides the equivalent contribution to the 
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ground state of two-electron ions. This follows from the fact that between s states, only the 
monopole part of the l/r 12 operator contributes. This radial contribution of the monopole 
part is exactly given by the potential in Eq. ( |B10| ). Moreover the retarded part of the 
Coulomb interaction vanishes in this case. Finally, the exchange correction only involves the 
spin of the two electrons and corresponds to a multiplication of the function F ls ls (Za) by 
two. 

The above arguments can be extended for all cases where an electron of arbitrary quan- 
tum number interacts with a s electron. The s electron couples only to the monopole term 
in the angular expansion of the Coulomb interaction. If, however, the electrons are not 
identical, then there is an exchange term with additional multipole terms. In this case, 
the retardation contribution to the Coulomb interaction is also non-zero. Obviously in the 
relativistic case the magnetic part of the electron-electron interaction should be considered. 

Because of these considerations we can compare our results to the Coulomb part of 
the calculation done by the Goteborg Group. In Fig. |^ we plot also the function F{Zoi) 
from Refs. P,|T5|]. The difference between the two calculations is displayed on Fig. The 



agreement is very good for medium-Z, while the difference between the two calculations 
increases with increasing Z, which is due to the inclusion of finite nuclear size in Refs. P,|T9" 



while the present results are for a point nucleus. This comparison thus provides the finite 
nuclear size effect on the two-electron self-energy. The difference at low Z (5 and 10) are 
due to numerical inaccuracies. Since no uncertainties are given in Ref. | 19[] , which contains 



more accurate values than Ref. [[|, we assume an uncertainty of 1 in the last digit (note 
that FiZa) in Ref. Jl9| is two times ours.) 



The low-Z behavior of the self-energy correction for to the Coulomb interaction to 



Fi s ,is(Za) is known from the work of Araki |22[] and Sucher |23" |, and 1/Z expansions from 
Drake IE31 to be 



F? s f s {Za) = \ 



2.588819 

^-21og2)log(2(Z«) 2 ) , (118) 
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while the magnetic part is 

Ful(Za) = I (119) 



We cannot directly compare our value for Fi Sj i s (Za) with those of Yerokhin et al. 
because their results also include the magnetic and retardation contribution to the self 
energy and our model does not. This correction contributes even at very low Z since it 
contains the free-electron anomalous magnetic moment from the vertex correction. The 
difference between the total contribution (for point nucleus) from Ref. Q and the present 



work is plotted on Fig. [LCJ for 20 > Z. Evidently, the magnetic interaction contribution to 
the two-electron self-energy is much larger than the contribution of the finite nuclear size. A 
simple fit of the difference between the present result and the one in Ref. j|] with a second- 
order polynomial yields 0.124 for the contribution of the anomalous magnetic moment, in 



good agreement with the value in Eq. ( |119|) . From the figure, it is evident that for Z as 



low as 5 the higher-order terms still make a significant contribution. On the same figure 



we also plot the Breit contribution from Ref. [19], which is in agreement with the difference 



between Yerokhin et al. and the present work for Z > 20, and matches reasonably well the 
extrapolated values even down to Z = 1. 

VII. CONCLUSION 

In this paper we describe a method of approximately evaluating two-electron radiative 
corrections that can easily be generalized to the direct evaluation of the correction repre- 
sented by the diagrams in Fig. 0. Accuracy and correctness of the method and programs 
is assessed by extensive comparisons with numerical derivatives of well-known one-electron 
self-energy results for a Coulomb perturbation. It is demonstrated that the method can work 
down to Z = 5 in some cases with reasonable accuracy. With the use of a more accurate 



Green's function evaluation and convergence acceleration techniques, following Refs. [20,21 



it is likely that calculation can be performed for He. The results presented in the present 
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paper also provides approximate self-energy screening corrections in any ion with less than 
10 electrons, thus providing a valuable, QED-based replacement for methods based on the 
Welton approximation |26| , |2T| or other, less efficient, screening schemes as used in atomic 
structure codes. It is also equivalent to the direct Coulomb contribution for some states of 
helium-like ions. 

This method could also be used with numerical Dirac-Fock potentials and wave functions 
from one of the codes in Refs. [p8p9f . Preliminary tests show that good numerical accu- 
racy can be achieved. Such an approach would provide more accurate self-energy screening 
corrections for the outer shells of very heavy transuranic elements or for inner hole binding 
energies 
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APPENDIX A: NUMERICAL TEST 

As a consistency check on the computer code, we carry out a test calculation in which 
the correction terms are generated by numerical differentiation of the unperturbed Coulomb 
self energy with respect to the nuclear charge Z . This should give the same result as the 
screening calculation where both the unperturbed potential and the perturbing potential 
are the Coulomb potential, with an appropriate normalization factor. 

In other words, we consider the potential 

(Z + AZ)a /k . 

V(x) = -- — (Al) 

Ob 

and let 
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V(x) = V (x) + SV(x) 



(A2) 



where 



V (x) 



Za 



x 



SV{x) 



AZa 



(A3) 



x 



If the level shift E(Z) is known as a function of Z for the Coulomb potential, then the exact 
correction due to 8V(x) is E(Z + AZ) — E(Z) and the first-order correction in AZ is 



Thus the first-order perturbation due to the potential 8V(x), with unit charge shift AZ = 1 
should be exactly equal to the derivative with respect to Z of the Coulomb level shift E'(Z). 

APPENDIX B: ORIGIN EXPANSION OF THE FIRST-ORDER CORRECTION 

TO THE WAVE FUNCTION 

This origin expansion is made with the use of the differential equation for the first-order 
correction to the wave function in Eq. (PH). However this expansion has a different form 
depending on whether one uses a Coulomb perturbing potential or a potential created by 
an other electron. 



In the case of a Coulomb perturbing potential, the correction to the wave function must 
have a logarithmic contribution, so as to cancel the a/x term in the lowest order of the 
development. We write 



SE(Z) = AZ lim 




AZ —E(Z) = AZE'(Z) 



(A4) 



1. Coulomb perturbation potential 



5fi(x) = ax' 



EC 0) Viog(*)£A, ( 



0)^1 

i ■ I 



(Bl) 
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and replace in Eq. (|91|), together with a series expansion of the unperturbed wave function, 
which behaves as uj near the origin. We then extract coefficients of x- 7 log (a;) and of x J and 
solve for the Q" and Aj coefficients. The coefficient of the log(x)/x term for an unperturbed 
wave function of angular symmetry K n is 

-7A? } + («« - A)A[, 2) = 



(B2) 



where the unperturbed wave function origin behavior is given by u = A — 1, A = Jk^ — 7 2 , 
7 = Z<x The determinant of this equation is zero and thus we can write 

A« = ^^ 2) . (B3) 

7 

The general equation for the term of order i is 

—K n — A - z)A? } + 7 Af } = (1 + K)AS 2 \ 

(B4) 

- 7 A, (1) + (« n - A - z)Af } = (1 - E n )\f\ 

The determinant of the linear system in Eq. ( |B4| ) is given by (2A + i)i and is nonzero for 
i > 0. By solving order after order, all higher-order terms can be expressed as a function of 

(2) 

Aq , 7 and K n . The expressions are all relatively simple since the unperturbed wave function 
does not have a logarithmic contribution. The non-logarithmic terms are obtained by first 
defining a series expansion for the unperturbed wave function 



°° I -\ T-' 



(B5) 



-j=o 3- 

where % = 1, 2 and iVj is a normalization factor. The equation derived from the term of 
order 1/x is given by 

-«„ - AJC, +7C0 - iV 2^2 H ^—\> 

(B6) 

-7Ci 1} + («n - A)d 2) = Arf - Aq 2) 



where we have used Eq. ( |B3| ) , and which again has a zero determinant. Explicit expressions 
of the n = 1 and n = 2 wave functions can be found in Ref. [|18|,|15|] (Note that in Ref. [FJ 
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the norm N3 of the 2p 3 / 2 wave function in Eq. (A3) should be iV 3 = 7 5 " 25 7[2r(5 - 25')}). 

(2) 

Requiring the compatibility of the two equations enables to calculate X as 

(2 ) _ K + A)iy iy 4 0) - Tivrf 



A, 



2A 

' an A A 1 ) 



(B7) 



One can then obtain a relation between Q and Q using one of the two equations in ( p6| 

„2a(i) 



so 



(1) 



7 2 Cr + K-A)Af- 7 iV 2 ^ 0) 



:bs) 



7(k„ + A) 

All ( coefficients can thus be expressed as a function of (q 2 \ These coefficients must be 
determined from the normalization condition of the perturbed wave function. This obliges 
to explicitly write n > 0, as a function of (q 2 ^ rather than keeping them as function of 
Cn-i- The latter expressions are simpler, but the former are very large. We use Mathematica 
to build the equations obeyed by the an \$ coefficients, evaluate the explicit expressions 
of extracting the part which depends on £g and the one which doesn't, and generating 
FORTRAN code. The code can have hundreds of lines for each piece of C?- The final 
expression of df can finally be recast as 



Sfi(x) 
5f 2 (x) 



ax 



ax 



, 0') xi 



s? ] (i + Ef =1 Cm) + ET=o C$ + log(x) Ef =0 AV : f 



-{j)x> 

j 



(J) 



v (i) 



(B9) 



A comparison for a small value of x of the expansion and of the value obtained by the use of 
the reduced Green's function as described in Sec. |111 G| yield two values of Q , one for each 
component of the wave function. A comparison of the two values provide a good check of the 
algebra. We compute Q? an d A,- coefficients up to j = 3. The value of obtained from 
each component of the wave function at x — 0.0005 agree with an accuracy of 13 significant 
figures, for n = 1 and n = 2, n n = —1, 1, —2. 



2. Electron screening potential 



The screening potential described in Eq. ( |C2| ) leads to an origin expansion rather different 
than the one described in the preceding section. In order to evaluate the origin expansion of 
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the first-order wave function correction we first evaluate the screening potential expansion. 
Using the origin expansion of the wave function and Eq. (|C2]), one can easily show that 

(BIO) 



SV (j \x) = 5V ij) (0) + x 2X ^ ]T 8V® 



\n=0 



X 



where Ay) = \/«y) — 7 2 (the origin behavior of the screening wave function is Ay) — 1) and 



8V® (0) =a dxx 
Jo 



fl{x) + f 2 2 {x) 



(Bll) 



The asymptotic expansion is obtained by substituting Eq. ( piO|) in the Poisson equation 
obeyed by the potential 

d 



dx 2 



(xSV®(x)) + ax fl{x) + f 2 (x 







(B12) 



and expanding the two radial component of the wave function in powers of x. 

With such an expansion the shape of the origin expansion of the first order correction to 
the wave function is 



8fi(x) = ax } 



X Jl 



j=0 

+x 2A o-)^aS j) - 



(B13) 



We obtain the equation for Aq by looking at the coefficients of x 2A ( J ' +A . We get for the 



2A m +A 



equation of order 1/x: 



,(°)jfT/(j) 



(-1 - 2A (i) - n n - A)A[, 1) + 7 Aj, 2) = -N 2 <p 2 "'6V ( 

-jx® + (-1 - 2\ w +K n - a)aJ 2) = -m^sv® 



(B14) 



-T rv n — /\j/\ — —xviyi uvq 

We note that in this case this equation in inhomogeneous and has a non-zero determinant. 
The equation for Q is obtained from the term of order 1/x 2 as 

(- Kn - \)<™ + 7 Cf 







7 Ci 1} + {K n ~ A)C 



(2) 



(B15) 
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APPENDIX C: MODEL POTENTIALS 



One of the models considered here for a screening potential is the spherically averaged 
potential that arises from the charge distribution of another electron in state j in the atom: 

1 



5V (j \x 2 ) 



a 



dQ 2 dxx- r 

47T J J \X 2 — Xl\ 



oo 



a I dxi ; 

/o max(i 2 ,a;i ; 



x 



fi,(j)( x i) + fl(.j)( x i) 



(CI) 



To facilitate numerical integration, this equation is written as 

*2 ~ / 1 1 



X 2 



a 

— ) — a 
x I JO 

x 



(C2) 



for x 2 < Xq, or as 



5V®(x 2 ) = — - a r dx x x\ (—- — 

Xo Jx2 \Xo Xi 



a-2 
X 



(C3) 



for x 2 > Xq, for a suitable value of xq. The expectation value in (IC4 ) is evaluated with the 
aid of the identity 

dE n 



dZ 



(C4) 



where E n is the energy eigenvalue of the screening wave function. 



The crossover point is taken to be xo —2/Jl — E^. The integral in ( |C2|) is evaluated by 
20 point Gauss Legendre quadrature with a new integration variable t over the range (0, 1) 
defined by X\ = x 2 t 4 , and the integral in (|C3|) is evaluated by 25 Gauss Laguerre quadrature 



with a new integration variable s over the range (0, 00) where X\ = x 2 +s/(2yl — E%). This 
prescription gives a precision of better than one part in 10 12 for the range 1 < Z < 100, as 
determined by comparing results of the two methods of integration in ( |C2"| ) and 
The corresponding derivatives are 
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-^8V<t\x 2 ) 
dx 2 



~ 2 / do&\$~^ 

Xn J 



(C5) 



for x 2 < Xq, or 



dx- 



-5V {j \ 



X 2 



a 

x 2 



X 1 



.'■2 



dx x x\ fl(j)(xi) +flu){xi) 



(C6) 



for X 2 > Xq. 

The derivatives are calculated with the same integration methods as the described above 
for the function 8V^\x). 

A simple additional model potential, useful for testing code, is generated by employing 
an exponential charge distribution, which corresponds to the replacement 



,0) 



(C7) 



in Eqs. (C2), flClf), (|C5|) , and ( P6|) and leads to the analytic potential 



SV cxp (x) = --a(- + r)i 
x \x 



—2-yx 



(a 



with 



«7 



(C9) 



and the derivative 



d_ 

dx 



5V exp (x) 



a 



x- 



+ a ( -r + — + 2 7 z ) e 



-27a; 



(CIO) 
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TABLES 

TABLE I. Comparison between the present direct calculation (Dir.) with 5V(r) = a/r and a 
calculation using numerical derivative of the one-electron self-energy (Num. Der.). Comparison is 
done for the scaled function F(Za) . 



Is 



z 


Num. Der. 


Dir. 


20 


-11.1383 (3) 


-11.1384 (8) 


50 


-6.35896 (5) 


-6.3589 (2) 


90 


-6.05030 (3) 


-6.0502 (5) 


2s 


Z 


Num. Der. 


Dir. 


20 


-12.2267 (8) 


-12.227 (1) 


50 


-8.0243 (1) 


-8.0243 (3) 


90 


-10.03340 (6) 


-10.0334 (6) 


2pi/ 2 


Z 


Num. Der. 


Dir. 


20 


0.3190 (8) 


0.319 (6) 


50 


-0.2444 (1) 


-0.2445 (8) 


90 


-2.30603 (4) 


-2.30607 (7) 


2p 3 /2 


Z 


Num. Der. 


Dir. 


20 


-0.606 (1) 


-0.61 (1) 


50 


-0.9056 (1) 


-0.906 (2) 


90 


-1.37474 (3) 


-1.3748 (3) 
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TABLE II. Self-energy screening function f(Za) for Is, 2s, 2pi/ 2 and 2p 3 / 2 electrons 



F(Za) for Is screened by 



Is 


2s 


2pi/ 2 


2p 3 /2 


5.171 (5) 


-0.491 (2) 


-0.243 (2) 


-0.242 (2) 


3.882 (1) 


-0.3681 (4) 


-0.1832 (5) 


-0.1809 (5) 


2.901 (1) 


-0.27670 (1) 


-0.14115 (1) 


-0.13555 (1) 


2.7386 (3) 


-0.26184 (10) 


-0.1345 (1) 


-0.1279 (1) 


2.4726 (3) 


-0.23772 (7) 


-0.12394 (8) 


-0.11536 (8) 


2.1726 (1) 


-0.21111 (4) 


-0.11286 (5) 


-0.10092 (5) 


2.0917 (1) 


-0.20410 (4) 


-0.11010 (5) 


-0.09696 (5) 


1.83507 (8) 


-0.18267 (3) 


-0.10245 (3) 


-0.08403 (3) 


1.73770 (6) 


-0.17506 (2) 


-0.10028 (3) 


-0.07891 (2) 


1.62045 (5) 


-0.16663 (2) 


-0.09870 (2) 


-0.07243 (2) 


1.55825 (4) 


-0.16272 (1) 


-0.09867 (2) 


-0.06875 (2) 


1.48502 (4) 


-0.15912 (1) 


-0.10011 (1) 


-0.06402 (1) 


1.43296 (4) 


-0.15809 (1) 


-0.10334 (1) 


-0.06003 (1) 


1.40906 (3) 


-0.158821 (9) 


-0.10661 (1) 


-0.057727 (10) 


1.3935 (2) 


-0.160737 (9) 


-0.11085 (1) 


-0.055657 (9) 


1.38620 (2) 


-0.166065 (8) 


-0.119409 (9) 


-0.052952 (8) 


1.39010 (7) 


-0.16998 (1) 


-0.124865 (10) 


-0.051754 (7) 


1.4214 (2) 


-0.18313 (4) 


-0.14160 (4) 


-0.04933 (2) 


1.43690 (2) 


-0.188136 (6) 


-0.147674 (8) 


-0.048726 (6) 



F(Za) for 2s screened by 

Is 2s 2pi /2 2p 3/2 

-7.951 (9) -3.162 (6) -4.180 (6) -4.171 (6) 

-5.770 (1) -2.2806 (8) -3.0141 (9) -2.9876 (9) 

-4.740 (1) -1.8624 (4) -2.4615 (5) -2.4120 (5) 
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-4.1691 (4) 


-1.6260 (2) 


-2.1503 (2) 


-2.0717 (2) 


-3.8535 (3) 


-1.4896 (1) 


-1.9720 (1) 


-1.8563 (1) 


-3.7137 (3) 


-1.4208 (1) 


-1.88388 (9) 


-1.71942 (8) 


-3.7207 (2) 


-1.4072 (2) 


-1.86978 (8) 


-1.63864 (9) 


-3.87690 (9) 


-1.44860 (6) 


-1.92989 (6) 


-1.60346 (6) 


-4.2184 (2) 


-1.55689 (8) 


-2.08125 (9) 


-1.61098 (8) 


-4.3155 (2) 


-1.58889 (8) 


-2.12572 (10) 


-1.61786 (8) 




F(Za) for 2p 1 / 2 screened by 




Is 


2s 


2P1/2 


2p 3 / 2 


0.10 (6) 


-0.03 (2) 


-0.06 (2) 


-0.06 (2) 


0.064 (5) 


-0.027 (2) 


-0.049 (2) 


-0.049 (2) 


-0.012 (1) 


-0.0373 (9) 


-0.0620 (6) 


-0.0611 (6) 


-0.119 (5) 


-0.056 (3) 


-0.088 (4) 


-0.084 (4) 


-0.2566 (5) 


-0.0817 (3) 


-0.1250 (2) 


-0.1158 (2) 


-0.4316 (2) 


-0.1164 (2) 


-0.1757 (1) 


-0.1547 (1) 


-u.uuu \ \y> j 


-0.1636 (7) 


-0.245 (7) 


fl 909 (1} 


-0.9751 (5) 


-0.2300 (6) 


-0.3446 (9) 


-0.2604 (5) 


-1.4367 (1) 


-0.33002 (7) 


-0.49508 (7) 


-0.33416 (6) 


-1.5564 (2) 


-0.3563 (2) 


-0.5348 (1) 


-0.35143 (9) 



Is 



F(Za) for 2p 3/2 
2s 



screened by 
2pi/2 



2p 3 /2 



-0.7 (2) 
-0.663 (3) 
-0.691 (4) 
-0.737 (1) 
-0.794 (1) 
-0.8599 (4) 



-0.21 (9) 
-0.185 (3) 
-0.189 (2) 
-0.1992 (3) 
-0.2154 (4) 
-0.2369 (2) 



-0.3 (1) 
-0.285 (3) 
-0.287 (2) 
-0.3004 (4) 
-0.3223 (4) 
-0.3524 (2) 



-0.3 (1) 
-0.281 (3) 
-0.279 (2) 
-0.2849 (4) 
-0.2962 (4) 
-0.3111 (2) 
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70 


-0.9330 (3) 


-0.2642 (2) 


-0.3912 (5) 


-0.3286 (3) 


80 


-1.0123 (3) 


-0.299 (1) 


-0.440 (1) 


-0.3486 (8) 


90 


-1.0972 (2) 


-0.3419 (2) 


-0.5022 (1) 


-0.3704 (1) 


92 


-1.114 (9) 


-0.351 (6) 


-0.516 (6) 


-0.374 (6) 



TABLE III. Comparison between F(Za) using the partial-wave renormalization method [|| 
and this work. More accurate numerical values from Ref. are used. i?rms is the nuclear 
mean-spherical charge radius used (in Fm). Since no error estimate is provided with Ref. fl9|l , in 
which 2F(Za) is tabulated we use an error of 1 on the last displayed figure. 



z 


i?rms 


Ref. [19| 


This work 


Diff. 


5 




-5.1745 (5) 


-5.171 (5) 


0.004 (5) 


10 




-3.8820 (5) 


-3.882 (1) 


0.000 (1) 


18 


3.423 


-2.89995 (5) 


-2.901 (1) 


-0.001 (1) 


24 


3.643 


-2.47240 (5) 


-2.4726 (3) 


-0.0002 (3) 


32 


4.07 


-2.09145 (5) 


-2.0917 (1) 


-0.0002 (1) 
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4.480 


-1.73720 (5) 


-1.73770 (6) 


-0.00050 (8) 


54 


4.78 


-1.55720 (5) 


-1.55825 (4) 


-0.00105 (7) 


66 


5.21 


-1.43035 (5) 


-1.43296 (4) 


-0.00261 (6) 


74 


5.37 


-1.38890 (5) 


-1.3935 (2) 


-0.0046 (3) 


83 


5.519 


-1.38080 (5) 


-1.39010 (7) 


-0.00930 (8) 


92 


5.860 


-1.41630 (5) 


-1.43690 (2) 


-0.02060 (5) 
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TABLE IV. Comparison of the present result for lithium-like ions with earlier work (eV). 
Results from Refs. [||-||,31-|34] all include finite nuclear size correction, not included in the present 
work. We use results without exchange. The difference with Ref. ||] on the 2p 1 / 2 screening is due 
to a programming error in the code used in 1991. 



Orbital screened by 



AE 



Ref. [§ DifT. Refs. |]Q Ref. Q Ref. § Ref. || 



Ref. [3 



Is 



2s 

2P1/2 
2P3/2 



Is 
2s 

2P1/2 
2P3/2 

Is 
Is 
Is 



-3.76607 (5) 
-0.49310 (2) 
-0.38705 (2) 
-0.12771 (2) 
-1.41385 (8) 
-0.50991 (5) 
-0.365 (3) 



-3.79 
-0.48 
-0.33 

-1.44 
-0.32 



0.02 
-0.01 
-0.06 

0.03 
-0.19 



-1.375 (5) -1.389 -1.39 -1.375 (30) -1.385 (15) 

-0.485 (5) -0.506 -0.505 -0.475 (30) -0.495 (15) 

-0.356 



Is 2 2s 
Is 2 2p 1/2 
Is 2 2p 3/2 



-11.3460 (1) -11.42 0.08 
-9.32606 (8) -8.88 -0.45 
-8.518 (4) 



Transitions (Valence+Core) 
2p 1/2 -* 2s 2.0200 (2) 

2p 3 /2 -* 2s 2.828 (4) 



2.55 -0.53 



2.01 (1) 



Transitions (Valence) 
2pi/ 2 -> 2s 1.8079 (1) 



1.78 (1) 



2P3/2 



2.098 (4) 



1.766 1.77 
2.065 



1.80 (6) 



1.78 (3) 
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FIGURES 




(A) (B) 

FIG. 1. Feynman diagrams for radiative corrections to the electron-electron interaction. (A) 
represents the wave function correction and (B) is the vertex correction. 

= ^ = r 




(A) (A') (B) 

FIG. 2. Feynman diagrams for radiative corrections to the electron-electron interaction after 

extraction of reducible (A) and irreducible (A') parts for the wave function correction (A) of Fig. 
1 





W SE 5V 

FIG. 3. Feynman diagrams for the self-energy perturbed by an external potential. In the case 

of the spherically-averaged screening potential used in this work, the diagram in the left originates 

in the diagram (A) in Fig. ||, and represents the wave function correction, while the diagram in the 

center comes from diagram (A') in Fig. |, using the (symbolic) relation ^G(E) = G(E) ■ G(E). 

The diagram on the right comes from the vertex correction. 
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Derivative 

Direct Calc. 

■ T. by T. 



-0.1496065 




10 15 10 13 10 11 10 9 10" 7 

Cut-off 

FIG. 4. Comparison between the convergence of the sum calculated using numerical derivation 
of the unperturbed self-energy sum as in Eq. (j80|) to the direct evaluation using perturbed wave 
function and reduced Green's function as defined in Eq. ( |101|) , and to an evaluation in which terms 
in Eqs. (|8l] ) to ( |84|) are derived for each n value in the sum (T.-by-t deriv.). Evaluation is done for 
a 2p 3/2 state and r = 0.9992, y = 1, t = 0.99, Z = 20. 
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25 



45 
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FIG. 5. Screening of a Is electron by a Is, 2s, 2pij% or 2p 3 / 2 electron. Sunnergren: Correction to 
the Coulomb part of the electron-electron interaction from Ref . , improved calculation following 
Ref. m. 
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10 30 50 70 90 



Z 

FIG. 6. Screening of a 2s electron by a Is, 2s, 2p±/ 2 or 2^3/2 electron. 
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FIG. 7. Screening of a 2pi/ 2 electron by a Is, 2s, 2pi/ 2 or 2p 3 / 2 electron. 
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FIG. 8. Screening of a 2p^/ 2 electron by a Is, 2s, 2pi/ 2 or 2p 3 / 2 electron. The large error bar 
for Z = 10 are only due to the low-energy part of the vertex correction 
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FIG. 9. Comparison between calculation using the partial- wave renormalization method H and 



this work. More accurate numerical values from Ref. [19| are used. 
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FIG. 10. Comparison between calculation for point nucleus (including Coulomb and magnetic 
part of the electron-electron interaction) from Ref. [|j and the Coulomb value from this work. The 
solid line represents a second-order polynomial fit to the difference and extrapolates to —0.124 for 
Z = 0, as expected from Eq. ( |119| ). Also plotted are the Breit values from Ref. [pTSj] , which are in 
reasonable agreement with the extrapolation. 
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